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Abstract 

Wing  rock  is  an  uncommanded  rolling  motion  which  has  been  observed  for  slender 
aircraft.  The  interaction  of  vortices  is  at  the  heart  of  the  phenomena.  Computational  fluid 
dynamics  can  be  used  to  predict  wing  rock  through  time  accurate  simulation.  However, 
a  parametric  search  based  on  unsteady  calculations  is  computationally  costly.  Different 
approaches  to  the  calculation  of  wing  rock  are  considered  in  this  reportr,  including  direct 
stability  calculations  to  determine  onset  angles,  linearised  time  domain  simulation  and 
reduced  order  modelling  using  proper  orthogonal  decomposition.  Finally,  the  formulation 
of  reduced  models  for  damping  and  LCO  prediction  is  given. 

This  work  has  been  funded  by  the  European  Office  of  Aerospace  Research  and  Devel¬ 
opment  under  contract  FA8655-03-1-3044.  This  report  is  the  final  report  for  phases  1  and 
2  of  this  contract. 


1  Introduction 

The  desire  for  increased  speed  and  agility  has  led  to  aircraft  designs  with  increasing  wing 
sweep  and  the  addition  of  highly  swept  leading  edge  extensions.  A  common  dynamic 
phenomenon  experienced  by  slender  aircraft  flying  at  high  angle  of  attack  is  known  as 
wing  rock.  A  classic  example  of  slender  wing  rock  occurred  on  the  Handley  Page  115 
research  aircraft  which  was  designed  to  study  the  aerodynamic  and  handling  qualities  of 
slender  aircraft  at  low  speed  [1].  At  angles  of  attack  above  20°  the  aircraft  first  experienced 
wing  rock,  with  the  maximum  roll  amplitude  experienced  being  around  40°  at  30°  angle 
of  attack.  This  motion  was  suppressed  by  reducing  the  angle  of  attack  or  applying  an 
input  to  the  aileron.  Aircraft  which  experience  such  motions  tend  to  have  highly  swept 
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surfaces  or  have  long  slender  forebodies  which  produce  vortical  flows.  At  some  critical 
angle  of  attack  the  aircraft  can  experience  a  roll  oscillation  which  grows  in  amplitude 
until  a  limit  cycle  oscillation  is  reached.  A  loss  in  roll  damping  is  usually  associated  with 
wing  rock.  Although  wing  rock  is  a  complicated  motion  which  involves  several  degrees  of 
freedom,  the  primary  motion  is  a  rolling  motion  around  the  aircraft  longitudinal  axis.  This 
motivates  experimental  and  computational  studies  involving  only  one  degree-of-freedom 
rolling  motions  to  adequately  reproduce  the  dynamics. 

Wing  rock  of  slender  delta  wings  has  been  studied  experimentally  over  the  past  two 
to  three  decades  and  several  review  papers  have  been  published  within  the  last  decade 
[2,  3,  4,  5].  It  should  also  be  noted  that  wing  rock  can  occur  for  nonslender  delta  wings 
[6]  and  even  rectangular  wings  of  low  aspect  ratio  (less  than  0.5)  [7]. 

Levin  and  Katz  [8]  performed  an  experimental  study  into  wing  rock  of  76°  and  80° 
sweep  delta  wings  with  the  wing  mounted  on  a  sting-balance.  It  was  observed  that 
only  the  80°  wing  would  undergo  self-induced  roll  oscillations  for  the  given  experimental 
conditions  (Reynolds  number,  bearing  friction  and  wing  moment  of  inertia),  therefore  it 
was  concluded  that  the  wing  aspect  ratio  must  be  less  than  1  for  wing  rock.  The  Reynolds 
number  based  on  the  root  chord  was  5  x  10'5  and  wing  rock  was  observed  for  an  angle  of 
attack  of  20°.  The  presence  of  vortex  breakdown  over  the  wing  was  found  to  limit  the 
amplitude  of  the  LCO.  During  the  free-to-roll  motion  a  loss  in  the  wing  average  lift  was 
observed  relative  to  the  static  lift  for  the  same  angle  of  attack.  It  was  also  observed  that 
side  forces  on  the  model  were  high,  indicating  that  the  presence  of  wing  rock  will  have 
a  strong  influence  on  free  flight  models.  Increasing  the  wind  tunnel  speed  increased  the 
amplitude  of  the  LCO  but  the  reduced  frequency  of  oscillations  remains  almost  unchanged. 

Katz  and  Levin  [9]  performed  an  experimental  study  of  a  delta  wing  /  canard  config¬ 
uration.  The  canard  has  an  aspect  ratio  of  0.7  with  the  wing  having  an  aspect  ratio  of 
1.  The  Reynolds  number  based  on  the  wing’s  root  chord  was  3  x  106  and  the  mechanical 
friction  was  minimal  (less  than  0.2  g-cm).  Unlike  the  results  of  Levin  and  Katz  [8]  where 
wing  rock  was  not  observed  for  this  identical  wing,  with  the  reduced  mechanical  friction 
in  the  experimental  setup  wing  rock  occurred.  The  effect  of  the  canard  was  observed  to 
increase  the  effective  leading  edge  sweep  of  the  configuration.  As  such,  with  the  increased 
effective  sweep,  the  75°  wing  /  canard  configuration  has  an  enlarged  wing  rock  envelope. 
Complex  vortex  interactions  were  present  (since  vortices  were  produced  by  the  canard 
and  the  wing)  and  nonsymmetric  oscillations  occurred. 

Arena  and  Nelson  [10]  undertook  a  series  of  comprehensive  studies  on  an  80°  sweep 
delta  wing  undergoing  wing  rock.  An  air  bearing  spindle  was  developed  for  the  free-to-roll 
tests  that  allowed  an  isolation  of  the  applied  torques  due  to  the  flowfield.  This  implies  that 
the  mechanical  friction  coefficient  is  effectively  zero.  Motion  history  plots  were  obtained, 
as  well  as  static  and  dynamic  flow  visualisation  of  vortex  position  and  breakdown  location, 
static  surface  flow  visualisation,  steady  and  unsteady  surface  pressure  distributions.  It 
was  observed  that  there  was  a  rate  dependent  hysteresis  in  vortex  location  due  to  a  time 
lag  in  the  motion.  This  time  lag  was  found  to  produce  a  dynamically  unstable  rolling 
moment  able  to  sustain  the  wing  rock  motion. 

Rolling  motion  about  the  longitudinal  axis  of  a  delta  wing  has  been  computed  using 
CFD  by  several  researchers  [11,  12,  13,  14].  Recently  a  comprehensive  numerical  study 
of  wing  rock  was  conducted  by  Saad  [15]  using  a  three  degree-of-freedom  flight  mechanics 
model  for  a  generic  fighter  aircraft  configuration  (forebody,  65°  leading  edge  sweep,  and 
vertical  fin).  Roll,  sideslip  and  vertical  degrees  of  freedom  were  allowed.  Including  the 
sideslip  degree  of  freedom  was  found  to  delay  onset  of  wing  rock  and  reduce  the  wing  rock 
amplitude. 

Time  domain  CFD  based  predictions  are  a  powerful  tool  for  studying  wing  rock.  How- 
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ever,  the  computational  cost  can  be  high.  A  way  of  reducing  the  cost  of  computing  para¬ 
metric  searches  for  aeroelastic  stability  behaviour  was  proposed  by  Morton  and  Beran 
from  the  US  Air  Force  Laboratories  [16].  Their  method  uses  dynamical  systems  theory  to 
characterise  the  nature  of  the  aeroelastic  instability,  with  this  additional  information  con¬ 
centrating  the  use  of  the  CFD.  In  this  way  the  problem  of  locating  a  one  parameter  Hopf 
bifurcation  was  reduced  from  multiple  time  marching  calculations  to  a  single  steady  state 
calculation  of  a  modified  system.  This  modified  system  calculates  the  value  of  the  param¬ 
eter  for  which  an  eigenvalue  of  the  system  Jacobian  matrix  crosses  the  imaginary  axis. 
A  convection-diffusion  problem  was  used  to  evaluate  the  approach  [17]  and  the  method 
was  then  applied  to  an  aeroelastic  system  consisting  of  an  aerofoil  moving  in  pitch  and 
plunge  [16].  The  linear  system  was  solved  using  a  direct  method  and  this  motivated  the 
use  of  an  approximate  Jacobian  matrix  to  reduce  the  cost.  Some  robustness  problems 
were  encountered  when  applying  the  method,  particularly  at  transonic  Mach  numbers.  A 
complex  variable  formulation  of  the  problem  was  introduced  in  [18]  which  resolved  some 
of  these  problems.  An  approach  considered  to  reduce  the  difficulties  of  applying  a  direct 
solver  to  large  linear  systems  was  to  use  domain  decomposition  to  reduce  the  size  of  the 
system  at  the  expense  of  an  outer  iteration  over  the  domains.  This  was  tested  on  the 
model  problem  in  references  [17]  and  [18]. 

The  problems  introduced  by  using  a  direct  solver  were  resolved  in  [19]  where  a  sparse 
matrix  formulation  was  used  to  make  feasible  the  solution  of  the  linear  system  for  much 
larger  grids.  The  Newton  iteration  was  modified  to  enhance  the  efficiency  of  the  scheme  fol¬ 
lowing  work  on  approximate  Jacobian  matrices  for  CFD  only  problems  [20].  The  method 
was  shown  to  be  effective  for  tracing  out  flutter  boundaries  for  symmetric  aerofoils  mov¬ 
ing  in  pitch  and  plunge,  with  reductions  of  two  orders  of  magnitude  in  the  computational 
time  required  when  compared  with  time  marching.  This  benefit  was  demonstrated  for  a 
three  dimensional  problem  in  reference  [21]. 

This  report  presents  the  formulation  of  both  a  nonlinear  and  linearised  time  domain 
method,  a  reduced  system  through  Proper  Orthogonal  Decomposition  (POD),  an  eigen¬ 
value  based  method  using  the  inverse  power  method  and  finally  a  direct  augmented  solver. 
A  comparison  is  made  of  these  methods  for  an  80  degree  delta  wing,  with  an  emphasis  on 
evaluating  the  computational  efficiency.  The  paper  continues  with  the  formulation  of  a 
solver  for  the  coupled  CFD-roll  equations,  followed  by  a  set  of  fast  methods,  direct  stabil¬ 
ity,  linearised  time  domain  and  proper  orthogonal  decomposition,  which  use  the  full  model 
to  predict  wing  rock  onset.  The  full  model  is  then  evaluated  for  a  test  case  involving  an 
80-degree  wing,  with  grid  dependency,  time  accuracy  and  modelling  influence  considered. 
Finally  the  fast  methods  are  compared.  The  formulation  of  reduced  models  based  on  the 
direct  solution  formulation  is  given  in  appendices. 


2  Formulation 

2.1  Flow  Solver 

All  simulations  described  in  this  paper  were  performed  using  the  University  of  Glasgow 
PMB  (Parallel  Multi-Block)  RANS  solver.  A  full  discussion  of  the  code  is  given  in  refer¬ 
ence  [20] .  PMB  uses  a  cell  centred  finite  volume  technique  to  solve  the  Euler  and  Reynolds 
Averaged  Navier-Stokes  (RANS)  equations.  The  diffusive  terms  are  discretised  using  a 
central  differencing  scheme,  and  the  convective  terms  are  discretised  using  Osher’s  ap¬ 
proximate  Riemann  solver  with  MUSCL  interpolation.  Steady  flow  calculations  proceed 
in  two  parts,  initially  running  an  explicit  scheme,  then  switching  to  an  implicit  scheme 
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to  obtain  quicker  convergence.  The  linear  system  arising  at  each  implicit  step  is  solved 
using  a  Krylov  subspace  method.  The  pre-conditioning  is  based  on  Block  Incomplete 
Lower-Upper  BILU(O)  factorisation  which  is  decoupled  across  blocks.  For  time-accurate 
simulations,  Jameson’s  pseudo-time  (dual-time  stepping)  formulation  [22]  is  applied,  with 
the  steady  state  solver  used  to  calculate  the  flow  steady  states  on  each  physical  time  step. 

For  the  RANS  simulation  the  two  equation  k -lu  turbulence  model  is  used  for  closure. 
It  is  well  known  that  most  linear  two-equation  turbulence  models  over-predict  the  eddy 
viscosity  within  vortex  cores,  thus  causing  too  much  diffusion  of  vorticity  [23].  This 
weakens  the  vortices  and  can  eliminate  secondary  separations,  especially  at  low  angles  of 
attack  where  the  vortices  are  weakest.  The  modification  suggested  by  Brandsma  et  al. 
[24]  for  the  production  of  turbulent  kinetic  energy  was  therefore  applied  to  the  standard 
k-u;  model  of  Wilcox  [25]  to  reduce  the  eddy-viscosity  in  vortex  cores  as 

Pk  =  min{P%,  (2.0  +  2.0min{0,  r  —  1  })p(3*kuj}.  (1) 

Here  is  the  unlimited  source  term  for  the  production  of  k,  and  r  is  the  ratio  of  the 
magnitude  of  the  rate-of-strain  and  vorticity  tensors.  When  turbulent  kinetic  energy  is 
over  predicted  in  the  vortex  core,  it  will  be  limited  to  a  value  relative  to  the  dissipation 
in  that  region.  This  modification  was  found  to  improve  predictions  compared  with  the 
standard  k-cu  turbulence  model  and  is  used  in  the  RANS  simulations  presented. 


2.2  Semi  Discrete  System  and  Jacobian 

The  coupled  CFD-rolling  motion  equations  are  written  in  semi-discrete  form  as 


(2) 


where 


and 


w  = 


W  a 
W/ 


R  = 


Ra 

% 


(3) 

(4) 


are  the  state  and  residual  vectors,  partitioned  into  aerodynamic  and  flight  variables.  Here 
wa  is  the  vector  of  CFD  solution  variables  in  conserved  variables  in  curvilinear  coordinates 
(i.e.  Vi  pi,  ViPiUi ,  ViPiVi ,  ViPiWi ,  UjE)  where  the  subscript  i  indicates  the  cell  number,  V)  is 
the  cell  volume  and  p,  u,  v,  w  and  E  denote  the  fluid  density,  three  components  of  velocity 
and  total  energy  respectively).  The  vector  Ra  is  the  CFD  residual  vector  resulting  from 
the  composition  of  the  cell  fluxes  calculated  using  Osher’s  method.  The  residual  in  each 
cell  depends  on  the  aerodynamic  variables  in  13  cells  (i.e.  that  cell  and  4  adjacent  cells  in 
each  curvilinear  direction).  The  residual  also  depends  on  the  grid  speeds  x*  and  the  pitch 
angle  of  the  model  a,  which  is  indicated  by  writing 


Rfl  =  R«(wa,Xi,a). 


(5) 


The  roll  equations  are 


^=R, 

dtt  1 


(6) 


where  wj  =  [cj)t,  (j)]T  and  R f  =  [/j6)c  +  D<pt,  4>t]T ■  Here  <f>  is  the  roll  angle  about  the  body 
axis,  Cic  is  the  rolling  moment  which  has  been  non-dimensionalised  as  Cjc  =  1/pooU^cf 
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where  and  Uoo  are  the  freestream  density  and  velocity  respectively,  cr  is  the  centreline 
chord,  D  is  the  coefficient  of  friction  and  p,  =  Poocr/Ixx  where  Ixx  is  the  rolling  moment 
of  inertia.  The  rolling  moment  depends  on  the  fluid  solution  but  for  a  fixed  fluid  solution 
(in  fact  for  a  fixed  fluid  pressure)  is  independent  of  cj>  and  ct.  Hence,  the  rolling  residual 
is  written  as 

R/  =  R/K'W/)  (?) 

with  the  dependence  on  ws  coming  only  from  the  second  component  cj).  The  coefficient 
of  friction  D  is  assumed  to  be  zero  for  the  physical  model  but  is  exploited  below  for  Hopf 
Bifurcation  calculations. 

The  full  equations  can  be  written  to  emphasise  the  dependencies  as 

H-o(wa,  4>i  (j>ti  c>0 
pCtc{  w0)  +  D(j)t  . 

fa 


d_ 

dt 


W  a 

fa 


The  stability  of  an  equilibrium  w  =  wo  which  satisfies  R(wo)  =  0  is  determined  by  the 
eigensystem  of  the  Jacobian  matrix  A  =  <9R/<9w.  The  matrix  can  be  written  out  in  block 
format  as 


-  &R a 

3R„ 

d  R„  1 

A  = 

<9w  a 

r1  dwa 

0 

d(j>t 

D 

1 

d<f> 

0 

0 

= 

Aaa 
_  Afa 

1  1 

C3  ^ 

(8) 


The  term  Aaa  is  a  large  sparse  matrix.  For  Afa  the  rolling  moment  depends  linearly  on 
the  surface  pressures,  which  in  turn  depend  linearly  on  the  values  in  the  first  two  internal 
cells  from  the  boundary.  The  term  Aaf  is  derived  from  noting  that  if  initially  (f>  =  0  and 
the  grid  is  at  xo,  and  the  body  axis  points  along  the  x-axis,  then  the  instantaneous  grid 
locations  are 


x  = 


1  0  0 

0  coscj)  — sinq i> 
0  sincj)  coscj) 


(9) 


The  grid  speeds  are  then 


x*  =  fa 


0  0  0 

0  —sincj)  —coscj) 
0  coscj)  —sincj) 


(10) 


The  terms  in  Aaj  are  then  calculated  through  finite  differences  by  incrementing  cj>  or  cj)t, 
evaluating  x  if  0  has  been  incremented  or  x*  if  fa  has  been  incremented,  updating  the 
body  boundary  conditions  and  then  recalculating  the  aerodynamic  residual  for  the  finite 
difference  calculation. 


2.3  Augmented  Solution  Based  on  the  Critical  Angle 

2.3.1  Augmented  System 

We  consider  the  stability  problem  when  the  angle  of  attack  a  is  considered  as  a  parameter. 
Note  that  the  angle  is  used  to  calculate  the  wing  location,  keeping  the  freestream  velocity 
vector  parallel  to  the  x-axis. We  are  interested  in  finding  the  onset  angle  for  the  wing  rock 
and  assume  that  stability  is  lost  through  a  Hopf  Bifurcation.  In  this  case  the  Jacobian 
matrix  A  has  a  pair  of  purely  imaginary  eigenvalues  at  the  critical  angle. 

At  this  angle  we  therefore  have 

A(Pi  +  iP2)  =  *o;(Pi  +  iP2)  (11) 
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where  Pi  and  P2  are  the  real  and  imaginary  parts  of  the  critical  eigenvector  and  iui  is 
the  purely  imaginary  eigenvalue.  The  eigenvector  is  normalised  against  a  real  vector  S  so 
that 

Sr(P1  +  iP2)  =  0  +  zl.  (12) 

At  the  bifurcation  value  for  a  we  therefore  have  a  system  of  conditions 

R 

APi  —  LUt*2 

R-aug  =  AP2+L0P1  =0.  (13) 

STPi 

stp2-i 

If  there  are  n  degrees  of  freedom  in  the  semi-discrete  system  state  vector  w,  then  equa¬ 
tion  (13)  has  3n+2  equations.  Equation  (13)  can  be  solved  for  the  augmented  vector  of 
unknowns 


where  w  is  the  equilibrium  at  the  bifurcation  point,  Pi  and  P2  are  the  real  and  imaginary 
parts  of  the  critical  eigenvector,  uj  is  the  frequency  of  the  instability  and  a  is  the  Hopf 
bifurcation  angle. 

The  system  of  equations  is  solved  using  Newton’s  method.  One  iteration  is  written  as 

'A  0  0  0 

AwPi  A  — Iuj  — P2  Aa  Pi 

AwP2  Iuj  A  Pi  AaP2  ~  <ug)  =  ~Kug  (15) 

0  ST  0  0  0 

0  0  ST  0  0 

There  are  three  problems  with  using  this  method.  These  are  now  dealt  with  in  turn. 

2.3.2  Forming  the  Augmented  Residual 

The  right  hand  side  vector  includes  the  terms  AP  \  and  AP2.  These  are  formed  from 
Residual  evaluations  through  a  matrix  free  evaluation 

=  R(w  +  eP)  -  R(w  -  eP)  (16) 

2e  '  ’ 

This  avoids  the  need  to  form  the  Jacobian  matrix  A  explicitly.  In  fact,  the  matrix  free 
product  is  used  for  the  contribution  of  the  term  Aaa. 

2.3.3  Forming  the  LHS  matrix 

Approximations  to  the  LHS  matrix  are  possible  as  long  as  the  iterations  converge  at  a 
reasonable  rate  to  the  correct  solution.  It  has  been  found  that  for  CFD  only  and  aeroelastic 
augmented  bifurcation  calculations  that  using  the  Jacobian  matrix  Aaa  corresponding 
to  the  first  order  spatial  scheme  leads  overall  to  a  more  efficient  scheme.  When  this 
approximation  is  made,  the  resulting  Jacobian  matrix  corresponding  to  equation  (8)  is 
denoted  by  A.  This  is  because,  although  quadratic  convergence  is  lost  and  so  more  quasi 


6 


Newton  steps  are  needed,  the  linear  system  is  much  easier  to  solve.  The  approximate 
Jacobian  involves  seven  non-zero  blocks  per  cell  in  the  grid.  The  second  Jacobian  Aw  is 
difficult  to  compute  efficiently  and  this  term  is  also  neglected  in  the  approximate  Jacobian 
driving  the  augmented  Newton  solver.  This  Jacobian  matrix  is  then  written  as 


dK, 


aug 


dvs 


aug 


AO  0 
0  A  —Iio 
0  Ioj  A 
0  ST  0 
0  0  ST 


0  0 

- P2  A*Pi 

Pi  AQP2 

0  0 

0  0 


(17) 


The  terms  AaPi  and  AaP2  are  formed  using  a  finite  difference  approximation.  The 
first  equation,  which  calculates  the  equilibrium  solution  at  the  current  angle,  is  decoupled 
from  the  critical  eigenvalue  equations. 


2.3.4  Forming  the  Exact  Jacobian  Matrix 

The  key  issue  is  how  the  convective  terms  in  the  Euler  equations  are  treated.  All  other 
terms  in  the  coupled  Jacobian  matrix  are  treated  either  analytically  or  by  finite  differences 
as  described  in  previous  reports.  We  will  therefore  focus  on  the  Jacobian  of  the  fluid 
residual  Ra  with  respect  to  the  fluid  unknowns  wa. 

The  residual  for  one  cell  in  the  grid  is  built  up  from  fluxes.  Following  the  usual 
approach  for  Riemann  solvers, 

£ i,j  =  wr) 

where  w;  =  wi(wi_2)j,  Wjj,  wi+1>j)  and  wr  =  wr(wi_ijJ-,  wi+ij,  wi+2j).  Here 

fjj  is  computed  using  Osher’s  [26]  approximate  Riemann  solver.  The  left  and  right  states 
are  computed  from  the  cell  values  using  MUSCL  interpolation.  For  the  cell  interface  we 
are  considering  there  are  four  contributions  to  the  Jacobian  matrix  arising  from 

<9wi_2j 

dki 

dwj-lj 

dkj 

dw  ij 

dki 

<9wi+i  ,j 

Now,  for  the  approximate  Jacobian  used  to  drive  convergence  of  the  Newton  iterations, 
the  following  approximations  were  made: 


dwj-2,j 

dfjj  _  d(jj 
dwi-ij  dw  1 
dfjj  _  dfjj 
3w  ij  dwr 


dvri+i,j 
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These  approximations  are  exact  for  a  first  order  spatial  discretisation  where  w /  =  w,_ij 
and  wr  =  Wjj.  The  calculation  of  the  terms 


dw  i 


and 


dwr 


is  non-trivial  but  has  been  coded,  tested  and  used  in  the  Glasgow  CFD  solver  since  1998 
[28].  These  are  exploited  to  calculate  the  exact  Jacobian  terms  for  the  second  order  spatial 
discretisation  by  using  the  Chain  rule 


dfjj  _  dfjj  dw i 
<9wj_2j  dw  I  3wj_2j 


dkj 

dw  ij 


dfjj  dw  i  dfjj  dw  r 

dw i  dw i_ij  dwr  dwi_i  j 

dfjj  dwi  dfjj  dwr 

dwi  dwj  j  dwr  dwi  j 


dfij  dfij  dwr 

dwl+X)3  dwr  dwi+1J 


The  new  coding  of  the  terms  arising  from  the  dependence  of  w;  and  wr  on  the  cell  values 
is  shown  in  appendix  A. 

Some  care  must  be  taken  at  boundaries  where  halo  cells  are  used  to  simplify  im¬ 
plementation.  The  halo  values  are  functions  of  the  internal  values  w&i  =  wji(wi,W2) 
and  Wfo2  =  wj,2(wi,  W2)  with  the  exact  relationship  depending  on  the  type  of  boundary. 
Applying  the  Chain  rule, 


dfi,  _  dfb  dwi  dib  dwr  d%  dwi  dwbi  dib  dwr  dwbi  dfi,  dwi  dwb2 

dw±  dwi  dw\  dwr  dw±  dwi  dwbi  dw\  dwr  dwbi  dw\  dwi  dwb2  dw\ 


dib  _  dib  dwi  dfb  dwr  dfb  dwL  dwbl  dfb  dwj  dwb2 

dw2  dwi  dw2  dwr  dw2  dwi  dwb\  dw2  dwi  dwb2  dw2 

The  dependence  of  the  halo  values  on  the  interior  values  leads  to  similar  extra  terms  from 
the  adjacent  interfaces  to  the  boundary  also. 

The  second  order  Jacobians  were  tested  by  forming  matrix  vector  products  against 
random  vectors  and  comparing  with  the  results  from  a  matrix  free  product. 


2.3.5  Initial  Values  for  the  Eigenvector 

The  convergence  of  the  Newton  iterates  requires  a  good  initial  guess.  We  assume  that 
we  have  an  estimate  for  a  and  to.  Then  a  steady  state  CFD  calculation  gives  wQ,  with 
ws  =  0.  Finally,  the  inverse  power  method  is  applied  to  the  approximate  Jacobian  A  or 
the  exact  Jacobian  A,  using  lo  as  a  shift  to  give  an  initial  guess  for  Pi  and  P2. 


2.3.6  Augmented  Solution  Based  on  the  Friction  Coefficient 

The  augmented  solver  described  above  suffers  from  a  potential  difficulty.  The  fluid  so¬ 
lution  is  strongly  dependent  on  the  angle  of  attack.  This  means  that  decoupling  the 
equilibrium  and  stability  equations  is  not  a  good  strategy.  A  structural  parameter  makes 
a  better  bifurcation  parameter  since  the  fluid  equations,  which  are  far  more  numerous, 
are  potentially  less  sensitive  to  a  structural  parameter.  We  therefore  consider  a  second 
approach  to  calculating  the  bifurcation  angle. 

Let  us  consider  the  case  where  the  angle  of  attack  is  fixed  and  the  structural  friction 
coefficient  D  is  regarded  as  the  bifurcation  parameter.  This  has  the  following  advantages 

•  The  fluid  and  structural  equilibrium  is  independent  of  D  and  therefore  the  first 
component  of  equation  (13)  decouples  from  the  other  equations.  Several  of  the 
approximations  made  in  obtaining  equation  (17)  are  therefore  not  necessary. 

•  The  terms  Aaa ,  Afa  and  Aaf  are  all  independent  of  D  and  so  can  be  calculated  once 
and  for  all  at  the  start  of  the  Newton  iterations. 

•  A  consequence  of  this  is  that  the  incomplete  L-U  decomposition  can  also  be  calcu¬ 
lated  at  the  start  of  the  Newton  iterations  and  then  held  fixed.  This  is  likely  to  be 
a  good  approximation  since  most  of  the  matrix  does  not  change  during  the  Newton 
iterations. 

•  Finally,  the  terms  AqP\  and  A0P2  are,  except  for  one  component,  equal  to  zero. 

The  onset  friction  coefficient  can  be  computed  through  the  following  scheme.  First,  the 
fluid  steady  state  is  computed.  Secondly,  the  inverse  power  method  is  applied  as  described 
above  to  provide  initial  values  for  the  real  and  imaginary  parts  of  the  eigenvector.  Thirdly, 
the  matrices  Aaa ,  Afa  and  Aaf  are  computed  and  fixed  for  the  rest  of  the  computation. 
The  augmented  Newton  iterations  are  then  driven  by 

A  — Ioj  —  P2  AjjP  | 
dRaug  _  IiO  A  Pi  ADP2 

<9w  aug  ST  0  0  0 

0  ST  0  0 

where  the  solution  variables  are 


(19) 


The  LU  decomposition  of  the  matrix  in  equation  (18)  is  calculated  at  the  first  iteration 
and  is  then  fixed  thereafter. 

Now,  the  model  friction  is  zero  for  a  free-flying  aircraft  and  so  basing  the  bifurcation 
calculation  on  the  onset  friction  is  not  very  useful  in  itself.  However,  this  calculation  can 
be  efficiently  made  at  a  number  of  angles  of  attack  and  the  value  of  onset  friction  which 
is  zero  computed.  When  this  happens  we  have  the  onset  angle  for  zero  model  friction. 

2.4  Linearised  Solver  and  POD 

A  simpler  version  of  the  time  domain  solver  which  is  likely  to  be  adequate  for  the  prediction 
of  stability  can  be  obtained  by  linearising  the  semi-discrete  system  about  a  steady  state. 
This  also  provides  a  starting  point  for  developing  a  reduced  order  model  using  proper 
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orthogonal  decomposition.  The  state  vector  w  is  written  as  w  =  wo  +  w '(t)  where  wo  is 
an  equilibrium  value  and  w \t)  is  the  perturbation  about  this  mean.  To  first  order 


dw' 

dt 


^  <9R  . 

R(wo)  +  —  W 


(20) 


and  since  wq  is  an  equilibrium,  R(wq)  =  0.  Hence, 


dw1 

dt 


Aw' 


(21) 


to  first  order. 

Now,  given  a  set  of  solutions  from  a  time  marching  calculation,  a  set  of  modes  can  be 
computed  by  proper  orthogonal  decomposition  (POD)  The  computed  modes  are  written 
in  matrix  form  as 

*  =  [<h\<h\ . l</>n]-  (22) 

The  projection  of  the  solution  onto  the  modes  is  w'  =  $w,  where  w  is  a  vector  of  modal 
coordinates.  By  using  the  projected  solution,  equation  (21)  can  be  written  as 


and  premultiplying  by  <f>r  gives 


_  dw  ^ 

<f> —  =  A<hw 
dt 


(23) 


Cl  w  rr 

— —  =  TTy4<hw  =  Aw.  (24) 

dt  v  ; 


The  matrix  A  is  the  modal  Jacobian  matrix  whose  components  can  be  formed  as 

tR(w0  +  e(j)j)  -  R(w0  -  e<j>j) 


Aij  — 


2e 


(25) 


Once  the  elements  of  A  have  been  computed,  requiring  2  residual  evaluations  for  each 
mode  retained,  the  stability  of  the  equilibrium  can  be  tested  for  the  a  used  to  train  the 
modes  using  the  small  dimension  linear  system  of  ODE’s  in  equation  (42). 

This  method  requires  the  calculation  of  mode  shapes  using  snapshots  from  time  march¬ 
ing  calculations.  The  formulation  is  based  on  primitive  variables  i.e.  the  mode  shapes  are 
calculated  in  terms  of  snapshots  expressed  in  terms  of  p,  u,  v,  w  and  p.  An  equation  for 
the  evolution  of  the  primitive  equations  can  be  derived  using  the  chain  rule  as 


dw  dp 
dp  dt 


R(w) 


(26) 


which  gives 


dp 

dt 


dp 

dw 


R(w)  =  R*(p). 


(27) 


The  solution  is  then  linearised  as  p  =  po  +  p’  (t)  and  then  p' 


equation  is 

^  =  $*TA$*p  =  A*  p 
dt  K  •  •  • 


<b*p.  Finally,  the  projected 


(28) 


where 


R*(p0  +  e</>*) 


-R*(po 

2e 


(29) 
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3  Optimisation  of  Linear  Solver 

3.1  Overview 

A  key  part  of  the  computational  work  for  the  methods  described  in  previous  reports  is 
the  solution  of  a  large  sparse  linear  system.  The  package  Aztec  [29]  was  used  for  this 
purpose.  Aztec  is  a  general  sparse  matrix  package  which  implements  Krylov  subspace 
solvers,  several  preconditioners  and  has  a  parallel  option.  A  general  block  structure  for 
the  matrix  can  be  accomodated. 

The  key  decisions  before  parallel  implementation  of  the  stability  prediction  methods 
concern  the  linear  solver.  In  particular,  the  preconditioner  used  is  a  block  incomplete 
LU  decomposition  (BILU)  [30]  which  is  sequential  in  nature.  Extra  approximations  are 
needed  to  obtain  a  method  which  is  effective  in  parallel.  Efforts  were  directed  first  at 
examining  the  sequential  performance  of  the  linear  solver.  Then,  a  parallel  version  was 
considered.  This  work  is  described  in  the  current  section. 

3.2  GCR  Algorithm 

Eisenstat,  Elman  and  Schultz  [27]  developed  a  generalized  Conjugate  Gradient  rnehtod 
that  depends  only  on  A  rather  than  ATA  and  is  called  the  Generalized  Conjugate  Resid¬ 
ual  (GCR)  algorithm.  Saad  and  Schultz  developed  the  Generalized  Minimal  Residual 
(GMRES)  algorithm  which  is  mathematically  equivalent  to  GCR  but  is  less  prone  to 
breakdown  for  certain  problems,  requires  less  storage  and  arithmetic  operations.  How¬ 
ever  GCR  remains  the  easier  algorithm  to  implement,  especially  in  parallel.  The  GCR 
algorithm  is: 


ro 

= 

b  — 

Ax  o 

Po 

= 

ro 

For 

o 

II 

1,2, 

. . .  ,  until  convergence. 

Do: 

a.j 

{rj,Apj) 

{ApjyApj) 

xj+ 1 

= 

Xj  +  ajPj 

r3+ 1 

= 

r3  ~  ajAPj 

flij 

= 

(Arj+1,APi)  f  ■  _  0 

<■ APi,Api )  ’  IOT  1  ~  U’ 

1,2, 

Pj+ 1 

= 

rj+ 1  +  Ei=0  PijPi 

Enddo 


(30) 


To  calculate  the  the  vector  Arj  and  the  previous  Apj' s  are  required.  The  number 
of  matrix  vector  products  per  step  can  be  reduced  to  one  if  Apj+\  is  calculated  by 

3 

Apj+ 1  =  Arj+ 1  +  Y  PijAPi  (31) 

i= o 

This  may  not  be  beneficial  if  A  is  sparse  and  j  is  large. 

A  restarted  version  called  GCR(m)  is  defined  so  that  when  the  iteration  reaches  step 
m  all  the  pj' s  and  Apj' s  are  thrown  away.  Since  GCR  and  GMRES  are  closely  linked 
it  might  be  possible  to  use  one  of  the  standard  GMRES  subbase  methods  to  improve 
convergence. 


3.3  BILU(k) 

For  the  Block  Incomplete  Lower  Upper  (BILU)  factorization  the  matrix  is  partitioned 
into  5x5  matrix  blocks  associated  with  each  cell  in  the  mesh.  The  use  of  this  blocking 
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reduces  the  memory  required  to  store  the  matrix  in  a  sparse  matrix  format.  For  clarity 
matrix  elements  refer  to  these  blocks  from  now  on. 

Consider  a  general  sparse  matrix  A  whose  elements  are  atJ ,  i .  j  =  1 , ,n.  A  general 
Incomplete  factorization  computes  a  sparse  lower  triangular  matrix  L  and  a  sparse  upper 
triangular  matrix  U  so  that  the  residual  matrix  R  =  LJJ  —  A  satisfies  certain  contraints, 
such  as  having  entries  in  a  prescribed  pattern.  A  common  constraint  consists  of  taking 
the  zero  pattern  of  the  L  JJ  factors  to  be  precisely  the  zero  pattern  of  A.  However  the 
accuracy  of  the  ILU(O)  incomplete  factorization  may  be  insufficient  to  provide  an  adequate 
rate  of  convergence. 

More  accurate  block  Incomplete  LU  factorisations  allowing  extra  terms  to  be  filled 
into  the  factorisation  are  often  more  efficient  as  well  as  more  robust.  Consider  Tipdating 
the  a,ij  element  in  full  Gaussian  Elimination  (GE).  The  inner  loop  contains  the  equation 

"  ij  —  CLij  "ik  "hj  ■  (32) 

If  1  evij  is  the  current  level  of  element  aVJ  then  the  new  level  is  defined  to  be 

levy  =  min(levjj,  levjfc  +  1  evkj  +  1).  (33) 

The  initial  level  of  fill-in  for  an  element  atJ  of  a  sparse  matrix  A  is  0  if  at]  ^  0  and 
oo  otherwise.  Each  time  the  element  is  modified  in  the  GE  process  its  level  of  fill-in  is 
updated  by  equation  33.  Observe  that  the  level  of  fill-in  of  an  element  will  never  increase 
during  the  elimination.  Thus  if  ^  0  in  the  original  matrix  A,  then  the  element  will 
have  a  level  of  fill-in  equal  to  zero  throughout  the  elimination  process.  The  above  gives 
a  systematic  algorithm  for  discarding  elements.  Hence  ILU(k)  contains  all  of  the  fill-in 
elements  whose  level  of  fill-in  does  not  exceed  k.  The  algorithm  is  given  by: 

For  all  non  zero  elements  aij  define  lev  =  0 

For  i  =  2, . . . ,  n  Do: 

For  j  =  1,2, ...  ,i  —  1  and  for  lev  a ij  <  k 
"ij  =  11  ij:  "jj 

"il  " il  "ij "j l  l  j  f  1,  .  .  .  ,TT  (34) 

Update  the  levels  of  fill  in  for  non  zero  a ij 

EndDo 

If  lev  atj  >  k  then  a ij  =  0 

EndDo 


3.4  Evaluation  of  Aztec 

The  starting  point  for  studying  the  linear  solver  is  the  Aztec  package  which  has  been  used 
to  generate  all  results  to  date.  As  a  test  case  we  use  the  system  Arx  =  b  where 


Ar 


A  Ioj 
—Ilo  A 


(35) 


Here  A  is  the  Jacobian  matrix  of  the  CFD  equations  plus  the  roll  equations.  This  system 
is  used  in  the  inverse  power  iterations  and  is  close  to  that  used  for  the  augmented  solver. 

Various  calculations  were  carried  out  using  a  test  matrix.  First,  the  value  of  u  was 
set  to  zero  to  obtain  a  system  with  just  A.  Secondly  the  problem  was  solved  with  a 
representative  value  for  oj  and  with  the  various  orderings  for  ILU  and  BILU  factorisations. 
Finally,  one  of  these  cases  was  rerun  with  a  second  order  Jacobian  matrix  for  A.  The 
results  are  summarised  in  figure  1. 

The  following  conclusions  were  drawn  for  the  solution  of  the  linear  system  with  Aztec: 
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•  The  augmented  linear  systems  are  significantly  more  difficult  to  solve  than  CFD  only 
systems. 

•  The  second  order  Jacobians  cannot  be  solved  with  zero  fill-in  in  preconditioners 
whereas  the  approximate  first  order  Jacobians  can. 


iteration 


Figure  1:  Convergence  histories  for  TFQMR  solution  of  augmented  system  using  several  pre¬ 
conditioning  options. 

To  allow  an  evaluation  of  the  performance  of  Aztec,  and  to  facilitate  easy  testing  of 
various  options,  an  implementation  of  a  Krylov  method  with  BILU  preconditioning  was 
written  in  MATLAB.  A  general  version  of  the  preconditioning  was  written  which  allows 
levels  of  fill-in. 

First,  tests  were  carried  out  using  u  =  0.  Figure  2  shows  the  convergence  using  different 
preconditioning  schemes  on  the  second  order  Jacobian  matrix.  Convergence  is  obtained, 
even  with  BILU(O)  in  contrast  to  Aztec.  BILU(l)  gives  a  rate  of  convergence  approxi¬ 
mately  equal  to  the  first  order  Jacobian.  Preconditioning  the  second  order  Jacobian  with 
a  preconditioner  from  the  first  order  Jacobian  works  well  and  saves  45%  in  storage. 

Figure  3  shows  the  convergence  of  BILU(O)  for  the  first  and  second  order  Jacobains  with 
non  zero  lu.  The  effect  of  adding  the  non-zero  lo  terms  is  minimal. 

The  conclusion  from  this  initial  investigation  was  that  the  Aztec  package  was  not 
performing  in  an  optimal  way  and  that  an  implementation  of  the  linear  solver  would  have 
to  be  coded  up.  Various  options  were  considered  and  these  are  described  in  the  remainder 
of  this  section. 
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Figure  2:  Convergence  histories  for  GCG  solution  of  augmented  system  with  uj  =  0  using 
several  preconditioning  options. 


Figure  3:  Convergence  histories  for  GCG  solution  of  augmented  system  with  u>  =  0  and  to  =  1.0 
using  BILU(O). 

3.5  Formulation 

The  test  problem  is  written  in  real  form  above.  An  alternative  is  to  solve  the  complex 
form  of  the  system,  written  as  ( A  —  iu)v  =  Acv  =  d.  The  matrix  A  is  further  partioned 
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as 


(36) 


A  = 


Aff  Afs 
A  sf  ^4*5 


using  the  notation  introduced  in  report  2.  The  real  and  complex  systems  then  have  the 
following  matrices 


Aff  Afs 

Asf  ASs 

—  iui 

O  ^ 

^  o 

(37) 

Aff 

Afs 

U)I 

0 

Asf 

—col 

0 

0 

Aff 

col 

Afs 

(38) 

0 

— ujI 

Asf 

A  number  of  different  preconditioning  options  for  the  real  and  complex  systems  were 
tested,  using  a  FORTRAN  implementation  (which  made  the  treatment  of  complex  vari¬ 
ables  easier). 

Three  different  methods  have  been  tested  from  a  Incomplete  factorization  of  all  of  the 
terms  in  the  matrix  to  a  ’’blocked”  matrix  were  a  large  amount  of  decoupling  has  been 
used. 


3.5.1  Method  1 

This  is  standard  BILU(fc)  on  the  unsimplified  real  Ar  or  complex  matrix  Ac. 


3.5.2  Method  2 


This  is  BILU(/c)  on  the  block  diagonal  of  either  Ac  or  Ar  i.e. 


Abr  = 

Aff 

0 

—  iuj 

'  i 

0  ’ 

0 

CO 

CO 
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I 

and 


Aff  0  0  0 

0  Ass  0  0 

0  0  Aff  o 

0  0  0  Ass 


(39) 


(40) 


3.5.3  Method  3 

This  extends  the  blocking  in  Method  2  by  also  including  the  blocking  of  the  multiblock 
grid.  This  means  Aff  loses  all  its  inter  block  connectivity  following  the  method  used  in 
the  Glasgow  CFD  only  solver. 

3.6  Results 

Table  1  shows  the  number  of  non  zero  5x5  blocks  required  for  the  different  methods  with 
a  first  order  Jacobian.  Each  complex  block  requires  twice  the  storage  of  a  real  block  due 
to  the  real  and  imaginary  parts.  Even  taking  this  into  account  the  complex  formulation 
always  uses  less  memory  than  the  real  one.  Also  the  scaling  of  the  memory  requirements 
is  increasingly  favourable  for  the  complex  formulation  as  the  level  of  fill-in  increases.  The 
real  formulation  of  method  2  results  in  a  singular  preconditioner  due  to  Ass  being  singular. 
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Table  2  shows  the  number  of  non  zero  5x5  blocks  required  for  the  different  methods 
with  a  second  order  Jacobian.  Comparing  them  with  the  first  order  Jacobian  with  level  3 
fill-in  in  the  second  order  Jacobian  requires  4  or  5  times  the  storage.  However  there  is  a 
much  larger  decrease  in  storage  requirements  as  terms  in  the  preconditioner  are  removed 
compared  to  the  first  order  case. 

Figure  4  shows  the  differences  between  using  different  size  bases  in  the  CGR  method 
for  a  first  order  Jacobian  and  a  level  zero  fill-in.  There  is  sometimes  a  plateau  when  the 
new  basis  is  started  due  to  the  simplistic  method  used  to  restart  -  just  drop  the  whole 
basis  and  start  again. 

Figure  5  shows  the  differences  between  using  different  levels  of  fill-in  with  the  restarted 
CGR  method  with  a  basis  of  60  and  a  first  order  Jacobian.  As  the  level  of  fill-in  increases 
the  number  of  iterations  to  convergence  drops.  This  is  mostly  due  to  the  reduction  in  the 
size  of  the  plateaus  after  the  restart  as  the  gradients  outside  of  these  regions  are  similar. 

Figure  6  shows  the  difference  between  methods  when  using  level  3  fill-in  and  60  bases 
in  the  CGR  solver  with  a  first  order  Jacobian.  Neglecting  Afs  and  Asf  has  very  little 
effect  on  the  overall  rate  of  convergence.  However  removing  all  the  terms  that  cross  the 
multiblock  mesh  boundaries  cases  a  large  drop  off  in  performance.  This  has  never  been 
seen  in  the  Aff  cases  and  since  only  two  equations  have  been  added  to  the  system  it  might 
be  be  possible  to  recover  the  Aff  convergence  rate. 

Figure  7  shows  the  differences  between  methods  with  a  basis  of  60  in  the  CGR  solver 
and  a  second  order  Jacobian.  For  the  second  order  Jacobian  level  2  fill-in  is  not  adequate 
for  convergence,  even  for  method  1.  Level  3  fill-in  allows  for  fast  convergence  for  method 
1  and  2  but  is  slower  for  method  3.  To  deal  with  this  we  will  consider  the  following.  If 
two  blocks  are  on  the  same  processor  then  that  information  will  be  retained.  This  means 
if  one  processor  is  used  the  preconditioner  is  defined  as  in  method  2  If  each  processor  just 
has  one  block  the  preconditioner  is  defined  as  that  in  method  3.  It  is  possible  that  the 
sharp  leading  edge  is  causing  the  problem  and  hence  careful  placement  of  the  blocks  using 
this  preconditioner  may  give  good  performance. 


Method  Real  or  Number  of  Non  Zeros  5x5  Blocks 


Complex  BILU(O)  BILU(l)  BILU(2)  BILU(3) 


1 

Real 

396518 

806558 

1590985 

N/A 

2 

Real 

N/A 

as  Preconditioner  is  singular 

1 

Complex 

175854 

304151 

511390 

902146 

2 

Complex 

151667 

278915 

485085 

874617 

3 

Complex 

141603 

247315 

402803 

689403 

Table  1:  Table  of  the  number  of  non  zero  in  the  preconditioner  for  the  first  order  Jacobian 


Method 

Real  or 
Complex 

Number  of  Non  Zeros  5x5  Blocks 
BILU(3) 

1 

Complex 

4273227 

2 

Complex 

4241749 

3 

Complex 

2741219 

Table  2:  Table  of  the  number  of  non  zero  in  the  preconditioner  for  the  Second  order  Jacobian 
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Figure  4:  Comparison  of  the  influence  of  size  of  the  Bases  used  in  the  GCR  method  -  first  order 
Jacobian,  method  1  preconditioning  and  no  fill-in. 


Figure  5:  Comparison  of  the  Real  and  Complex  formulations  for  Method  1  using  a  first  order 
Jacobian. 
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Figure  6:  Comparison  of  the  influence  of  the  preconditioning  method  using  the  Complex  for¬ 
mulation  and  three  levels  of  fill-in  for  a  first  order  Jacobian  matrix. 


Figure  7:  Comparison  of  the  influence  of  the  preconditioning  method  using  the  Complex  for¬ 
mulation  with  a  second  order  Jacobian 
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3.7  Conclusions  on  Linear  Solver 

The  following  conclusions  have  been  drawn 

•  The  complex  form  is  more  reliable  and  uses  less  storage  and  hence  will  be  used 

•  The  facility  to  change  the  fill-in  on  the  fluid  Jacobian  is  useful 

•  We  can  now  solve  the  system  using  the  Jacobian  from  the  second  order  spatial 
scheme,  the  requirement  for  which  was  discussed  in  the  introduction 

•  Method  3  provides  a  simple  way  of  implementing  the  preconditioner  in  parallel 

•  The  FORTRAN  implementation  is  fast  compared  with  the  general  AZTEC  code  or 
the  Real  C  implementation. 


4  Parallel  Implementation 

4.1  Data  Decomposition 

The  data  decomposition  is  done  by  storing  whole  blocks  of  the  multiblock  grid  on  just 
a  single  processor.  This  means  partitioning  of  the  grid  can  be  consided  as  a  problem  of 
partitioning  the  blocks.  Each  processor  stores  a  certain  number  of  blocks  in  the  grid  and 
their  associated  fluid  variables.  The  structural  variables  are  treated  in  a  different  way. 
Each  processor  stores  all  the  structural  information.  This  is  currently  not  too  expensive 
as  the  number  of  structural  equations  is  small  compared  to  the  number  of  fluid  equations, 
but  might  have  to  be  revisited  for  different  problems. 

4.2  Decomposition  of  the  residual  and  the  Jacobian 

The  residual  and  Jacobian  are  decomposed  in  the  method  shown  in  figure  8.  Here  4 
processors  are  assumed  and  the  schematic  in  the  figure  indicates 
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Processors  with  no  solid  boundary  conditions  will  have  zeros  in  the  Asf  part  of  the  Jaco¬ 
bian. 

4.3  Matrix  vector  and  Matrix  Transpose  vector  operations 

Since  the  Jacobian  matrix  is  calculated  in  4  distinct  parts  Aff,  Afs,  Asf  and  Ass,  the 
matrix  vector  product  is  done  in  4  parts  which  are  summed  together.  The  products  with 
Afs  and  Ass  vector  are  trivial  as  the  structural  unknowns  are  on  all  processors.  The 
product  with  Asf  requires  a  “reduce  all”  at  the  end  of  the  multiply  to  get  the  information 
to  all  the  processors. 

A  different  ordering  is  used  for  the  fluid  Jacobian  of  Aff.  The  cells  are  partitioned  as 
shown  in  figure  9,  with  3  different  sets  of  cells.  The  internal  set  can  be  updated  without 
any  communication  and  hence  it  is  possible  to  complete  this  part  of  the  product  while 
waiting  on  the  messages  to  arrive  from  other  processors.  The  border  set  are  the  cells 
owned  by  this  processor  but  which  require  some  communication  with  other  processors  to 
be  fully  calculated.  The  external  set  is  not  updated  at  all  on  this  processor  but  is  used  to 
update  the  border  points.  Since  the  renumbering  of  Aff  has  modified  both  the  row  and 
column  indices  the  multiplying  vector  must  be  reordered  to  match.  The  output  product  is 
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Figure  8:  Effective  data  decomposition  of  the  Jacobian  and  right  hand  side  vector 

also  in  the  reordered  state  and  must  be  permutated  back  before  doing  the  preconditioner 
operations. 

The  reason  why  this  reordering  is  only  done  to  the  fluid  Jacobian  is  because  only  the 
matrix  vector  product  that  requires  communication  to  be  able  to  complete.  The  above 
renumbering  allows  for  the  maximum  amount  of  work  to  be  allocated  between  sending 
the  halo  messages  and  having  to  use  the  halo  information  from  another  processor. 

As  can  be  seen  in  figure  10  the  matrix  transpose  vector  operation  -  needed  to  compute 
the  left  eigenvector  of  the  system  -  looks  very  similar  to  the  matrix  vector  operations. 
The  only  difference  is  that  for  A the  processors  now  have  complete  columns  of  the  fluid 
variables  and  not  rows.  The  Ajj  product  requires  no  communication  to  form  the  vectors 
but  the  external  set  must  be  communicated  at  the  end  to  complete  the  operation. 

4.4  Load  Balancing 

To  maximise  the  parallel  preformance  of  the  code  the  blocks  must  be  carefully  assigned  to 
the  processors.  If  there  is  a  large  number  of  blocks  it  is  not  difficult  to  get  approximately 
the  same  workload  onto  each  processor  for  the  CFD  only  case.  However  in  the  IPM  the 
workload  for  non  solid  wall  fluid  points  is  greater  than  that  for  any  other  point  -  due  to 
the  contributions  from  Asf.  Hence  the  number  of  fluid  points  on  each  processor  and  the 
number  of  solid  wall  points  on  each  processor  should  be  balanced.  Lastly  the  number  of 
processor  block  boundaries  needs  to  be  kept  as  small  as  possible.  Not  only  does  this  keep 
the  size  of  the  messages  smaller  but  much  more  importantly  BILU(k)  can  be  calculated 
with  no  extra  appoximations  on  block  boundaries  if  both  blocks  are  on  the  same  processor. 

4.5  Preconditioner 

Method  3  is  used  for  implementing  the  preconditioner  in  parallel.  The  preconditioner 
contains  no  inter  block  boundary  information  and  hence  each  block  does  not  rely  on  any 
infomation  from  any  other  block.  If  this  preconditioning  does  not  provide  convergence  it 
its  possible  to  fill-in  the  block-block  interface  information  if  both  blocks  are  on  the  same 
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processor.  This  however  does  mean  the  convergence  of  the  linear  system  would  depend 
on  the  partitioning  of  the  blocks  to  the  processors. 


4.6  Scaling  Results 

The  results  in  table  3  were  obtained  by  solving  a  linear  system  arising  from  a  second 
order  Jacobian.  The  speedup  on  2  processors  is  very  good  considering  the  small  size  of 
the  problem.  The  speedup  on  4  processors  is  lower  due  to  a  small  5%  load  inbalance  and 
a  large  number  of  faces  needed  to  be  communicated  to  obtain  this  load  balance.  These 
results  are  therefore  considered  very  satifactory  and  open  up  the  capability  to  tackle  fine 
mesh  densities  using  parallel  processing. 


Processors 

Time  in  Seconds 

Speed  up 

1 

821 

2 

419 

1.95 

4 

237 

3.46 

Table  3:  Table  of  speedups  for  IPM  on  a  problem  with  22,000  cells 


5  Test  Case  and  Time  Marching  Results 

5.1  Test  Case  Description 

There  is  a  considerable  amount  of  published  experimental  data  for  80°  sweep  delta  wings. 
Hence  this  sweep  angle  was  selected  for  the  current  study.  The  geometry  of  the  wing  was 
identical  to  that  used  by  Arena  and  Nelson  [10].  The  wing  has  a  flat  upper  and  lower 
surface  with  a  45°  windward  bevel  and  a  root  chord  of  0.4222m.  The  moment  of  inertia 
for  this  wing  was  given  as  Ixx  =  0.00125  Kg-m2.  The  experiment  was  performed  at  a 
Reynolds  number  of  1.5  x  10°.  Since  PMB  is  a  compressible  flow  solver,  a  Mach  number 
of  0.2  was  used  to  avoid  convergence  difficulties  at  the  very  low  Mach  numbers  used  in 
the  experiments.  The  freestream  air  density  is  used  to  non-dimensionalise  the  moment 
of  inertia  of  the  wing  for  the  computations.  The  freestream  air  density  in  experiment  is 
unavailable  and  so  the  value  was  assumed  to  be  1.23 Kgm~3  (sea- level  ISA  conditions). 

5.2  Steady  State  CFD  Tests 

Before  simulating  wing  rock,  the  accuracy  of  steady  solutions  was  verified  with  a  grid 
refinement  study.  A  fine  grid  for  Euler  simulations  was  created  with  approximately  1.6 
million  points.  From  this  grid  two  levels  were  extracted  by  removing  every  second  grid 
point  in  each  direction.  The  RANS  solution  was  performed  on  a  grid  with  a  high  resolution 
in  the  vortical  and  boundary  layer  regions.  This  grid  has  1.8  million  points.  Steady 
state  solutions  were  computed  for  each  grid  with  the  residual  being  reduced  6  orders  of 
magnitude.  The  test  case  chosen  for  validation  purposes  has  the  wing  at  30°  angle  of 
attack  and  rolled  +10°. 

The  upper  surface  pressure  distributions  from  all  three  Euler  grid  levels  are  shown  in 
figure  11  for  the  chordwise  station  at  60%cr,  and  are  compared  with  the  measurements 
of  Arena  [31].  It  should  be  noted  that  the  experimental  results  indicated  laminar  flow  on 
the  upper  surface  of  the  wing  which  has  the  effect  of  moving  the  primary  vortices  inboard 
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Figure  11:  Grid  refinement  study  -  Surface  pressure  distributions  at  60%cr.  Lines  indicate  the 
coarse,  medium  and  fine  Euler  results. 

and  upwards  off  the  surface  of  the  wing.  This  should  be  kept  in  mind  when  considering 
wing  rock  which  is  attributed  to  hysteretic  vortex  movement  [10].  Examination  of  the 
surface  pressure  distributions  from  each  grid  indicates  that  the  solutions  are  not  grid 
converged.  However  Euler  simulations  of  delta  wing  flows  are  known  to  be  highly  sensitive 
to  grid  density  due  to  the  mechanisms  for  the  generation  of  vorticity  which  are  based 
on  the  truncation  error  of  the  discretisation,  [32]  and  therefore  the  current  results  are 
as  expected.  Despite  the  requirement  to  chose  an  appropriate  level  of  grid  density,  it 
is  known  that  Euler  simulations  can  realistically  predict  the  dynamic  response  of  leading 
edge  vortices.  Comparing  the  surface  pressure  distributions  at  60%cr ,  there  is  a  reasonably 
good  agreement  between  the  solution  from  the  medium  grid  and  experiment.  The  effect 
of  viscosity  can  be  seen  as  a  drop  in  the  primary  vortex  suction  peaks,  which  is  due  to 
the  presence  of  a  secondary  separation  shifting  the  vortex  location.  As  such  the  RANS 
suction  peaks  and  locations  compare  well  with  experiment. 


5.3  Wing  Rock  Response 

The  increase  in  pressure  difference  between  the  port  and  starboard  sides  as  the  grid  is 
refined  will  induce  a  stronger  restoring  moment  which  is  likely  to  produce  a  stronger  wing 
rock  response.  This  is  seen  in  figure  12  where  clearly  the  largest  amplitude  wing  rock 
is  predicted  by  the  fine  grid  (followed  by  the  medium  and  coarse  grids).  The  wing  rock 
amplitudes  for  the  fine,  medium,  and  coarse  grids  are  70°,  50°,  and  15°,  with  the  wing 
rock  amplitude  observed  in  experiment  being  40°. 

A  comparison  of  the  rolling  moments  from  the  medium  grid  and  experiment  through 
one  complete  steady  wing  rock  cycle  is  shown  in  figure  13.  Despite  the  maximum  roll  angle 
being  exceeded  in  the  Euler  solution  it  is  clear  that  there  is  reasonably  good  agreement 
with  experiment,  with  the  rolling  moment  distribution  and  magnitudes  being  predicted 
well.  As  discussed  by  Arena  and  Nelson  [10]  there  is  a  clockwise  loop  where  energy  is 
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Figure  12:  Grid  refinement  study  -  Wing  rock  histories 

added  to  the  system,  and  two  anti-clockwise  damping  lobes  when  energy  is  dissipated. 
Comparing  the  width  of  the  loops  where  energy  is  added,  more  energy  is  being  added  to 
the  system  in  the  Euler  simulations  in  comparison  to  experiment.  The  damping  lobes  are 
also  larger  in  the  Euler  solution  in  comparison  to  experiment. 

In  order  to  verify  the  temporal  accuracy  of  the  solutions,  a  time  step  refinement  study 
was  conducted.  It  was  found  that  the  non-dimensional  time  step  of  0.1875  provides  an 
adequate  temporal  resolution  for  the  current  wing  rock  studies.  This  equates  to  around 
120  time  steps  per  cycle.  For  forced  motion  studies  as  few  as  50  time  steps  per  cycle 
are  sufficient  for  temporal  convergence,  however  since  in  free-to-roll  studies  errors  amplify 
with  time  (where  auto-rotation  was  observed  to  occur  when  the  time  step  was  too  large), 
a  smaller  time  step  is  required  to  accurately  predict  the  wing  rock  behaviour. 

Roll  angle  histories  for  Euler  and  RANS  simulations  of  the  80°  sweep  delta  wing  at 
30°  angle  of  attack  are  shown  in  figure  14.  Comparing  first  the  wing  rock  amplitudes,  the 
amplitudes  from  the  RANS  solution,  Euler  solution,  and  experiment  are  35°,  50°,  and  40° 
respectively.  Clearly  the  effect  of  viscosity  is  to  reduce  the  wing  rock  amplitude.  This 
is  as  expected  based  on  the  previous  discussion.  With  the  higher  more  outboard  suction 
peaks  from  the  Euler  solutions  there  is  a  much  larger  wing  rock  response  in  comparison 
to  the  RANS  solutions,  where  the  suction  peaks  are  lower  and  more  inboard.  Comparing 
the  period  of  the  wing  rock  cycle  we  can  see  that  the  Euler  and  RANS  solutions  produce 
similar  results.  As  such  it  is  concluded  that  the  Euler  solutions  predict  qualitatively  the 
correct  vortex  dynamics  and  therefore  a  realistic  wing  rock  response. 

5.4  Onset  Angle 

Simulations  of  free-to-roll  motion  at  15°  and  20°  angle  of  attack  were  conducted  on  the 
medium  grid  using  the  Euler  equations.  The  roll  angle  histories  for  both  simulations  are 
shown  in  figure  15.  At  15°  angle  of  attack  the  solution  is  dynamically  stable  (i.e.  the  am¬ 
plitude  of  the  oscillations  decreases  due  to  aerodynamic  damping).  However  at  20°  angle 
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Figure  13:  Comparison  of  rolling  moments  through  a  steady  wing  rock  cycle  -  CFD  and  exper¬ 
iment 


Figure  14:  Effect  of  modelling  -  Comparison  between  RANS  and  Euler  solutions  for  a  =  30° 
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Figure  15:  Wing  rock  roll  angle  histones  -  a  =  15°  and  20° 

of  attack,  the  initial  roll  angle  of  10°  initiates  a  wing  rock  response  with  the  amplitude  of 
the  motion  increasing  with  time.  The  onset  of  wing  rock  observed  in  experiment  was  22° 
[10]. 

6  Fast  Methods  for  Predicting  Onset 

6.1  Overview 

The  simulations  described  above  involve  the  time  domain  simulation  of  the  wing  response 
to  an  initial  roll  perturbation.  Multiple  time  domain  calculations  to  determine  the  on¬ 
set  angle  for  the  instability  are  potentially  costly.  It  is  therefore  of  interest  to  consider 
methods  that  retain  the  full  modelling  fidelity  of  the  CFD  but  are  computationally  more 
efficient.  Various  methods  were  described  earlier  in  the  paper;  Hopf  Bifurcation  calcula¬ 
tion,  inverse  power  method,  linearised  solver  and  proper  orthogonal  decomposition  (POD) 
reduced  order  modelling.  In  this  section  these  methods  are  applied  to  the  wing  rock  pre¬ 
diction  and  their  efficiency  is  evaluated.  This  is  done  for  the  coarse  Euler  grid,  using 
first  and  second  order  spatial  discretisations.  The  results  of  the  previous  section  suggest 
that  qualitatively  the  behaviour  on  this  grid  is  correct,  and  that  methods  which  can  be 
generalised  to  solving  the  RANS  equations  on  fine  grids  will  be  of  use  for  quantitative 
predictions. 

6.2  First  Order  Spatial  Results 

It  is  useful  to  start  with  a  first  order  spatial  discretisation,  since  for  this  case  A  =  A. 

The  full  nonlinear  and  spatially  linearised  systems  should  be  able  to  predict  the  on¬ 
set  angle  through  the  growth  or  decay  of  the  response  to  an  initial  perturbation.  The 
responses  were  calculated  with  a  reduced  time  step  of  0.125  at  three  different  angles  of 
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Figure  16:  Convergence  of  onset  angle  for  the  direct  calculation  (first  order). 

attack,  namely  20,  24  and  27  degrees.  Results  for  the  linearised  and  nonlinear  methods 
are  identical.  The  response  of  the  roll  angle  is  damped  at  20  degrees,  grows  at  27  degrees 
and  is  near  to  being  neutrally  stable  at  24  degrees.  By  looking  to  the  ratio  of  consective 
peaks  in  the  traces  it  is  possible  to  estimate  that  the  onset  angle  is  24.3  degrees. 

The  inverse  power  method  can  be  used  to  calculate  the  behaviour  of  the  critical  eigen¬ 
value  as  a  function  of  the  angle  of  incidence  for  the  first  order  spatial  scheme.  A  shift  of  0.2 
was  found  to  result  in  convergence  to  the  critical  eigenvector.  The  eigenvalue  was  found 
to  converge  to  six  decimal  places  within  7  steps.  At  24  degrees  the  critical  eigenvalue  is 
—  1.65  x  10-4  +  z0.237  and  at  25  degrees  it  was  3.46  x  —4  +  -10.247.  Linear  interpolation 
suggests  that  the  eigenvalue  crosses  the  imaginary  axis  at  24.3  degrees. 

The  direct  calculation  based  on  the  onset  angle  was  given  an  initial  angle  of  20  degrees. 
From  this  the  first  order  spatial  problem  converged  rapidly  in  20  steps,  as  shown  in 
figure  16.  The  converged  onset  angle  is  24.3  degrees.  The  converged  frequency  is  0.239, 
which  is  also  in  agreement  with  the  inverse  power  method  and  the  linearised  time  domain 
predictions. 

The  direct  calculation  based  on  the  damping  coefficient  was  run  at  angles  of  21,  23 
and  25  degrees.  The  convergence  of  the  onset  damping  at  these  angles  is  shown  in  figure 
17.  Linear  interpolation  between  these  values  suggests  that  the  damping  goes  to  zero  at 
an  angle  of  24.3  degrees. 

Finally,  snapshots  were  collected  during  one  cycle  of  the  linearised  time  marching 
calculations.  From  these  POD  modes  were  generated  using  the  matlab  script  given  in 
appendix  A.  The  linearised  system  was  then  projected  onto  these  modes  to  provide  a 
small  order  linear  system.  Six  modes  were  retained  for  this.  The  eigenspectrum  of  the 
reduced  Jacobian  at  each  angle  was  calculated  in  matlab.  Again,  this  indicates  that  at 
20  degrees  the  system  is  stable,  at  24  degrees  it  is  just  below  neutrally  stable  and  at  27 
degrees  unstable.  By  linear  interpolation  of  the  real  part  of  the  critical  eigenvalue  it  is 
predicted  that  the  onset  angle  is  24.4  degrees  and  the  frequency  is  0.237. 


6.3  Second  Order  Spatial  Results 

Next,  the  second  order  spatial  scheme  was  used. 
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Figure  17:  Convergence  of  onset  damping  for  the  direct  calculation  (first  order). 

Time  marching  results  are  shown  in  figure  18.  Again  the  nonlinear  and  linearised 
schemes  give  identical  results.  The  response  at  23  degrees  grows  whereas  the  one  at  21 
degrees  decays.  Linear  interpolation  of  the  amplification  gives  an  onset  angle  of  22.4 
degrees. 

The  inverse  power  method  was  used  to  calculate  the  wing  rock  onset  angle.  3  Levels 
of  fill-in  were  used  for  the  linear  solver  which  allowed  the  residual  of  the  linear  system  to 
be  driven  down  11  orders  in  less  than  60  iterations. 

At  22  degrees,  after  7  inverse  power  iterations  the  critical  eigenvalue  had  converged 
to  6  significant  figures.  The  value  is  —1.381245  x  10-4  ±  iO. 223497.  At  23  degrees  the 
converged  value  is  1.805895  x  10~4  ±  iO. 232095.  These  angles  bracket  the  onset  angle, 
which  can  be  estimated  by  linear  interpolation  as  22.4  degrees. 

The  direct  calculation  to  find  the  onset  angle  did  not  convergence.  This  appeared 
to  be  due  to  the  sequencing  between  the  eigenvalue  and  equilibrium  equations.  The 
problem  was  not  present  for  the  first  order  spatial  system,  possibly  because  in  that  case 
the  eigenvalue  equations  were  hit  with  a  Newton  step  each  time,  which  resulted  in  a  high 
level  of  convergence.  The  direct  calculation  based  on  damping  was  made  at  angles  of  21 
and  23  degrees.  The  convergence  of  the  damping  in  each  case  is  shown  in  figure  19.  By 
linearly  interpolating  between  the  damping  values  predicted  the  onset  angle  is  estimated 
to  be  22.4  degrees. 

Snapshots  were  again  collected  during  one  cycle  of  the  linearised  time  marching  calcu¬ 
lations  at  21  and  23  degrees.  The  reduced  Jacobian  has  a  critical  eigenvalue  of  —0.0004  + 
70.218  and  0.00018  +  z0.232  at  21  and  23  degrees  respectively. 

The  cost  of  the  second  order  calculations  is  summarised  in  table  4.  The  linear  and 
nonlinear  time  marching  results  are  based  on  computing  five  cycles  of  the  response.  The 
reduced  model  was  based  on  modes  computed  from  one  cycle  of  the  linearised  solver. 
The  inverse  power  solver  is  relatively  more  costly  than  the  bifurcation  solver  due  to  the 
increased  cost  of  solving  the  full  as  opposed  to  the  approximate  linear  system. 


Figure  18:  Response  for  linear  and  nonlinear  time  marching  (second  order). 


Figure  19:  Convergence  of  onset  damping  for  the  direct  calculation  with  a  second  order  spatial 
scheme  (second  order). 


Method 

Relative  CPU  time 

steady  state 

1.00 

Bifurcation  (Damping) 

1.36 

Inverse  Power 

6.00 

Linearised  Time  domain 

13.10 

Nonlinear  Time  Domain 

81.00 

POD  reduced  model 

3.60 

Table  4:  Calculation  times  (second  order)  normalised  by  the  steady  state  calculation  time. 
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7  Conclusions 


A  range  of  methods  for  predicting  wing  rock  onset  have  been  presented.  The  basis  for 
the  predictions  is  a  residual  calculation  for  the  Euler  equations.  An  assessment  of  the 
influence  of  viscous  effects,  grid  dependency  and  time  accuracy  was  first  made.  The  fast 
methods  considered  are  direct  bifurcation  calculation,  linearised  time  domain  solution  and 
reduced  order  modelling  based  on  proper  orthogonal  decomposition. 

The  following  conclusions  have  been  drawn 

•  The  methods  provide  consistent  results  for  the  onset  angle 

•  The  bifurcation  calculation  for  the  onset  angle  was  unsatisfactory  due  to  sequencing 
problems  between  the  mean  flow  and  eigensystem  equations. 

•  An  alternative  approach  based  on  structural  damping  proved  successful  and  allowed 
the  onset  angle  to  be  calculated  efficiently. 

•  The  relative  cost  of  the  various  methods  showed  the  efficiency  of  the  direct  calcula¬ 
tion, with  the  linearised  time  domain  solver  being  an  order  of  magnitude  faster  than 
the  nonlinear  time  domain  solver. 
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A  Formulation  of  Reduced  Order  Model 

The  prediction  of  bifurcation  points  is  only  one  part  of  an  analysis  which  is  likely  to  be 
required.  The  pre  and  post  bifurcation  behaviour  is  probably  of  much  more  interest.  For 
example,  flight  tests  monitor  damping  characteristics  and  for  a  theory  to  be  of  use  in 
reducing  risk  and  cost,  the  pre-bifurcation  damping  must  be  readily  available.  There  is 
also  much  interest  in  limit  cycle  behaviour  post  bifurcation. 

This  implies  that  a  fast  method  to  predict  the  bifurcation  point  is  not  enough.  In 
addition,  methods  are  required  to  predict  behaviour  for  parameter  values  in  the  proximity 
of  the  bifurcation  point.  If  these  methods  revert  to  the  full  order  system  of  equations  then 
much  of  the  benefit  of  the  fast  onset  prediction  will  be  lost.  However,  the  structure  of 
the  dynamical  system  close  to  the  bifurcation  point  can  be  exploited  to  produce  small 
order  models  which  can  predict  the  important  dynamics  and  which  are  computationally 
efficient. 

In  this  section  the  formulation  of  two  approaches  to  be  examined  is  given,  along  with 
background  information. 

A.l  Centre  Manifold  Reduction 

Consider  the  nonlinear  system  of  equations 

w  =  R(w,a),  weT,  a  €  (41) 

where  R  is  sufficiently  smooth.  We  assume  that  we  are  at  a  Hopf  bifurcation  at  a  =  0 
and  hence  the  Jacobian  matrix  SR/<9w  has  2  and  only  2  critical  eigenvalues  with  zero 
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real  part  and  the  remaining  ns  =  n  —  2  eigenvalues  have  negative  real  parts.  Then  the 
system  (41)  can  be  transformed  to 


=  Ac(a)wc  +  g(wc,ws,a) 
=  ^4s(a)ws  +  h(wc,  ws,  a) 


(42) 


where  wc  6  5ft2  and  ws  G  5ftnc.  7lc(0)  is  a  2  x  2  matrix  with  its  eigenvalues  on  the 
imaginary  axis  and  7^(0)  is  a  nc  x  nc  matrix  with  no  eigenvalues  on  the  imaginary  axis. 
The  functions  g  and  h  have  at  least  quadratic  terms.  The  system  42  has  a  family  of 
invariant  manifolds,  locally  representable  for  small  ||a||  as 


wa  =  {(wc,  ha{ wc,  a))  :  wc  £  Jf2 


Wr 


<e} 


(43) 


where  e  >  0  is  sufficiently  small  and  ha  :  5ft2  x  5ftTO  — >  5ftns  is  smooth.  Moreover,  ha( wc,  0)  = 
hc( wc),  i.e.  Wq  coinsides  with  the  centre  manifold  of  Wc  the  system  with  no  paramater 
variation. 

The  Reduction  Principle  says  system  (42)  is  locally  topologically  equivalent  near  the 
origin  to 

=  Ac(a)wc  + g{wc,ha{wc,a),a) 

=  7ls(0)ws 

The  important  thing  to  notice  is  that  the  equations  for  wc  and  ws  are  decoupled  in 
equation  (44).  The  first  equation  is  the  restriction  of  equation  (42)  to  its  centre  manifold. 
The  dynamics  of  the  structurally  unstable  system  (42)  are  essentially  determined  by  this 
restriction,  since  the  second  equation  in  (44)  is  linear.  For  a  Hopf  bifurcation  with  (Ag2  = 
T'icu)  then  the  system  looks  like 


w  r 


Wo 


(44) 


/  o  -u  W  wci  \  (  G\(wci,wC2,  ws,  a) 

y  W  o  J  \  Wc2  J  y  G2(wcl,wC2,ws,a) 

=  Asws  +  H(wci,wC2,ws,a) 


(45) 


Due  to  the  tangent  property  of  Wq  we  can  parameterise  for  small  |a|  via  a  means  of 
a  projection  of  onto  Tc  =  {ws  =  0}.  For  a  first  approximation  we  shall  assume  a  has 
no  effect  on  this  parameterization. 

It  is  possible  to  rewrite  this  in  complex  form  by  use  of  the  variable  z  =  wc\  +  iwc 2  to 
obtain 

J  z  =  iuz  +  G(z,  z,  ws) 

W,  =  Asws  +  H(z,z,  ws) 


(46) 


where  G  and  H  are  smooth  complex-valued  functions  of  z,  z  G  C 1  of  at  least  quadratic 
order.  The  centre  manifold  Wc  can  be  locally  represented  as  a  graph  of  a  smooth  function 


Wc  =  {(z,  ws)  :  ws  =  h(z)} 

h  maps  5ft2  — >■  5ftn^2  and  due  to  the  tangent  property  of  Wc,  h(z )  =  0||wc||2.  The  Centre 
manifold  Wc  therefore  has  the  representation 

1  1 

w s  =  h(z,  z)  +  -W20Z2  +  Wnzz  +  -W02 Z2  +  0(|z|3),  (47) 

with  the  coefficients  Wij  G  C 2 .  Since  ws  must  be  real,  w\\  is  real  and  W20  =  wq2-  Using 
Taylor  expansions  in  z,  z,  and  ws  the  system  (46)  can  be  rewritten  as 

!z  =  iuiz  +  \G2QZ2  +  G\\zz  +  \G  02Z2 

+  \G2\z2z  +  (G10,  ^s)z  +  (G01,  ws)z  +  . . .  (48) 

ws  =  ^lsws  +  \H2qz2  +  Huzz  +  \Hq2z2  +  ... 
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where  G20,  Gn,  G0 2, 
and  H2 o  =  #02- 


Gioj 


G01,j 


G-2i  €  C1  and  Goi,  Gun  Hjj  €  Cn  2.  Since  w.s  is  real  Hu 
Qj+k 


G]k  dzidzk 


G(z,  z,  0) 


s  j  +  k  >  2, 


z=0 


d2 


dwSjdz 


G{z,z,  ws) 


,  j  =  1,2,  ...,n-2, 


2  — OjWg  —0 


g- 


dw,jdlG(z’*’w’) 


,  j  =  1,2,  ...,n-2, 


z=0,ws=0 


Qj+k 


H:,k  dzidzk 


^H(z,z,0) 


,  j  +  k  —  2, 


z=0 


is  real, 

(49) 

(50) 

(51) 

(52) 


On  substituting  the  representation  of  the  centre  manifold  in  (48)  and  equating  co- 
efhents, 

w2o  =  (2 iul  -  Asy1H2o 

wu  =  -C-'Hu  (53) 

=  (-2 iul  -  As)-lH2Q 


Where  /  is  the  identity  matrix  and  the  matrices  (2iujl  —  z4s),  As  and  (—2 iuil  —  As )  are 
invertible  since  0  and  d=2io;  are  not  eigenvalues  of  As. 


A. 2  Projection  Based  Reduction 

The  centre  manifold  theory  in  principle  can  be  used  to  reduce  the  system  to  two  degrees  of 
freedom  in  the  vicinity  of  the  bifurcation  point.  However,  there  are  considerable  difficulties 
of  implementation  when  using  systems  of  large  degree  and  functions  of  the  complexity 
encountered  in  the  wing  rock  project.  To  overcome  this  a  method  is  proposed  whose  main 
objective  is  to  maintain  the  original  structure  of  the  problem  as  much  as  possible,  allowing 
Jacobian  matrices  of  the  original  system  to  be  used  as  much  as  possible. 

The  full  system  can  be  transformed  by  using  only  the  vectors  corresponding  to  the 
critical  eigenvalues  of  A  and  its  transpose  AT .  These  are  calculated  by  using  the  direct 
method  for  the  onset  angle.  The  system  is  projected  onto  its  critical  eigenspace  and 
complement.  Suppose  we  have  a  Taylor  expansion  of  the  residual  function  R  about  the 
equilibrium  solution  wo  and  parameter  at  the  bifurcation  point  (Iq,  giving 

w  =  4w  +  F(w,/i),  w  €  (54) 

where  F(w ,  /t)  has  at  least  quadratic  terms  and  w  =  w  —  wo,  //  =  //  —  /Mj-  The  matrix  A 
has  a  pair  of  complex  eigenvalues  on  the  imaginary  axis  Ai  ,2  =  iuj,  lu  >  0.  Let  q  be  the 
right  eigenvector  corresponding  to  Ai.  Then  q  is  the  right  eigenvector  corresponding  to 
X2  and 

Aq  =  icoq,  Hq  =  —  icaq 
The  left  eigenvector  p  has  the  same  property 

A7  p  =  —iojp,  Ar p  =  iojp. 

These  can  be  normalised  such  that  (p,  q)  =  1  where  (p,q)  =  X)?=i  PiQi-  The  eigenspace 
S  corresponding  to  is  two  dimensional  and  is  spanned  by  {JJq,  9q}-  The  eigenspace 
T  corresponds  to  all  the  other  eigenvalues  of  A  and  is  n  —  2  dimensional.  Then  y  G  T  if 
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and  only  if  (p,  y)  =  0.  Since  y  G  while  p  is  complex  then  two  real  constraints  on  y 
exist  and  hence  it  is  possible  to  decompose  any  wGPas 

w  =  zq  +  5q  +  y 

where  z  G  C1,  zq  +  zq  £  S,  and  y  G  T.  The  complex  variable  z  is  a  coordinate  of  S  so 

I  z  =  (P,w) 

\y  =  w  -  (p,  w)q  -  (p,w)q 
since  (p,  q)  =  0.  So  the  equation  (54)  has  the  form 
f  i  =  ivz+(p,F(zq  +  zq  +  y,fa)) 

[  y  =  Ay  +  F(z q  +  zq  +  y,fi)-  (p,  F(zq  +  zq  +  y,  p)) q  -  (p,  F(zq  +  zq  +  y,  p))q 
This  system  is  (n  +  2)  dimensional  but  we  have  two  constraints  on  y. 


A. 2.1  Simplified  Damping  Model 

Now,  in  general  at  this  stage  a  centre  manifold  reduction  should  be  used  to  obtain  a 
relationship  between  y  and  z  which  allows  the  critical  dynamics  to  be  calculated  from  the 
z  equation  only.  This  treatment  allows  nonlinear  features  such  as  Limit  Cycle  Oscillations 
to  be  calculated  from  the  reduced  model.  A  description  of  the  method  to  perform  this 
reduction  is  given  in  appendix  A,  and  is  referred  to  in  the  results  section  as  the  centre 
manifold  reduced  model.  Here  we  are  interested  in  calculating  damping  for  parameter 
values  below  the  bifurcation  point.  In  this  case  the  influence  of  the  component  y  from  the 
non-critical  space  is  damped  faster  than  the  critical  component  z.  We  therefore  neglect 
the  influence  of  y  altogether  which  removes  the  need  for  the  central  manifold  reduction. 
Further  justification  for  this  approximation  will  be  given  from  results  for  a  test  problem 
below. 

Now,  in  general  at  this  stage  a  centre  manifold  reduction  should  be  used  to  obtain  a 
relationship  between  y  and  z  which  allows  the  critical  dynamics  to  be  calculated  from  the 
z  equation  only.  This  treatment  allows  nonlinear  features  such  as  Limit  Cycle  Oscillations 
to  be  calculated  from  the  reduced  model.  A  description  of  the  method  to  perform  this 
reduction  is  given  in  appendix  A,  and  is  referred  to  in  the  results  section  as  the  centre 
manifold  reduced  model.  Here  we  are  interested  in  calculating  damping  for  parameter 
values  below  the  bifurcation  point.  In  this  case  the  influence  of  the  component  y  from  the 
non-critical  space  is  damped  faster  than  the  critical  component  z.  We  therefore  neglect 
the  influence  of  y  altogether  which  removes  the  need  for  the  central  manifold  reduction. 
Further  justification  for  this  approximation  will  be  given  from  results  for  a  test  problem 
below. 

The  damping  is  therefore  determined  by  solving  the  equation 

i  =  iwz  +  (p,  F(zq  +  5q,  ft))  . 


This  system  is  two  dimensional.  Finally,  we  need  to  calculate  the  form  of  F.  Expanding 
the  function  R  in  a  Taylor  series  about  the  equilibrium  solution  wo  and  parameter  /jq , 
noting  that  in  fact  the  bifurcation  parameter  is  the  angle  of  attack  a,  gives 


<9R 

R(w,  a)  =  R(wo,  ao)  +  Aw  +  L>(w,  w)  +  C(w,  w,  w)  +  — —  a 

da 


1  <92R 


2  da 2 


d  T  AqA  T  —  Aaad  w  -)- 


1  dB(  w,  w) 


da 


-a. 


(55) 

(56) 
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A. 2. 2  Full  Reduction 


The  full  reduction  requires  the  application  of  the  centre  manifold  theory  to  bring  in  the 
influence  of  the  components  y  on  the  critical  coordinates  in  the  two-degree-of-freedom 
model. 

This  system  is  Taylor  expanded  in  z,  z  and  y  to  give  the  following  approximation 


\  z  —  iujz  +  \G20z2  +  G\\zz  +  \Gq2%2  +  \G2\z2z  +  (Gun  y)z  +  (Gm,  y)z  +  . . . 

\  y  =  Ay+  \h20z2  +  Huzz  +  H02z 2  +  . . . 

(57) 

where  G20,  G11,  Gq2,  G2\  €  C1;  G01,  G10,  Hij  €  Cn.  This  approach  leads  to  the 
scalars  and  vectors  all  being  functions  of  F  or  inner  products  of  p  and  F  (i.e.  the  original 
function).  The  manipulation  of  this  system  is  feasible,  even  for  systems  of  large  dimension. 
The  centre  manifold  can  be  represented  by 

y  =  h(z,  z)  =  ljh20z2  +  hnzz  +  \h02z2  +  O |z|3 


with  the  constraint  (p,  ht])  =  0.  The  vectors  hij  €  Cn  can  be  found  from  the  linear 
equations 

f  (2  iiol  —  A)h-2Q  =  H2q 

<  —  Ahn  =  H\\  (58) 

[  (—2 iul  —  A)h()2  =  Hq2 

These  equations  are  invertable  since  0,  and  ±2 iio  are  not  eigenvalues  of  A.  We  can  now 
write  the  restricted  equation  as 

I  =  hoz  +  \G2qz2  +  Guzz  +  \Gq2z2 

+  5(^21  —  2(Gio,  A  lH\\)  +  (G01 ,  (2iluI  —  A)  1H2q))z2z  +  . . . 

If  we  write  F(w)  in  terms  of  functions  B{x,y)  and  C(x,y,z)  defined  as 


then 


B(wi,w2) 


1  d2R 

2  dw2 


W1W2 


C(wi,w2,w3) 


1  d3R 

7W1W2W3 

o  aw 5 


F(w)  =  B(  w,  w)  +  C(w,  w,  w)  +  0||w||4 


Then  we  can  express 


(Gi  o,y)  =  (p,B(q,y)),  (G0i,y)  =  (p,R(q,y)) 

and  hence  the  restricted  equation  is  in  the  form 

I  =  iuiz  +  \G2qz2  +  Guzz  +  ^G02z2 

+  \{G2 1  -  2(p,  B{ q,  A^Hu))  +  (p,  B{ q,  (2 iujl  -  A)~1H2q)))z2z  +  . . . 

where 


G20  =  (p,  B{ q,  q))  Gn  =  (p,  B{ q,  q))  G02  =  (p,  B( q,  q))  G2i  =  (p,  C(q,  q,  q)) 


H20  =  B( q,  q)  -  (p,  B{ q,  q))q  -  (p,  B( q,  q))q 
Hu  =  B( q,  q)  -  (p,  B{ q,  q))q  -  (p,  B( q,  q))q 
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Using  these  equations  and  the  identities 


A  xq  =  — q  A  xq  =  — q  (2 iul  —  A)  xq 
iu  iu 

the  restricted  equation  can  be  rewritten  as 


—  q  (2  iul  -A)  Xq  =  -|-q 

IU  SlU 


Z  =  iuz  +  ~ff20^  +  911ZZ  +  ]^go2Z2  +  \g2lZ2Z 

where 

520  =  (p,  #(q,  q))  5n  =  (p,  q)) 

and 


521  =  (P,  0(q,  q,  q)) 

-  2(p,  B( q,  A~1B( q,  q)))  +  (p,  B( q,  (2 iul  -  A)~1B{ q,  q))) 

+  ^(p,5(q,q))(p,JB(q,q)) 

-  —  l(p,-B(q,q))l2  -  ^-|(p,^(q,q)}|2 

IU  6%U 

This  leads  to  the  following  algorithm: 

1.  calculate  /i,  P,  Q  and  u 

2.  calculate  i?(q,  q)  and  i?(q,  q) 

3.  solve 

f  (2 iul  -  A)h2o  =  B( q,  q) 

I  Ahu  =  B{  q,q) 


4.  form  C(q,  q,  q) 

5.  form  520,  5n,  521 

6.  time  march  ^  ^  ^ 

k  =  iuz  +  -g2oz2  +  5ii2:z  +  -g02Z2  +  -g2iz2z 


(59) 


A. 3  Perturbation  Method 

The  perturbation  method  of  multiple  scales  can  be  applied  to  determine  behaviour  close 
to  a  bifurcation.  This  approach  was  used  in  [35]  for  an  two  dimensional  aeroelastic  prob¬ 
lem  and  was  presented  in  [36]  for  the  Duffing  oscillator  and  in  [37]  for  a  three  equation 
model  problem.  The  idea  is  use  a  perturbation  parameter  e  to  separate  out  events  at 
different  time  scales.  The  terms  of  0(e)  provide  a  description  of  the  linear  dynamics  of 
the  system,  and  higher  order  terms  provide  information  about  the  effect  of  nonlinearities. 
The  presentation  of  the  method  follows  that  given  in  [35]. 

Write 

w  =  w0  +  ewi(T0,  T2)  +  e2w2(T0,  T2)  + .  (60) 

where  Tq  =  t.  and  T2  =  e2t.  Then,  expanding  the  residual  in  a  Taylor  series  about  the 
bifurcating  equilibrium  wo,  and  equating  terms  in  e  gives  the  following  relations 

Wl  =  @peiujT°  +  Ope~iu}T° 

w2  =  r  +  r  (6i) 
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0  =  A/310  +  4/?2©20 


where 


T  =  2zo0©  +  2z2©2e2ia;To 

Pi  =  qTAMp 

P2  =  qT  (2B(jp,  z0)  +  B(p,  z2)  +  ^(p.p.p)) 
Az0  =  ~7}B(  p,p) 


The  operator  A,t  denotes 


(2 iul  -  A)z2  =  -B( p,p). 


d2R 


(62) 

(63) 

(64) 

(65) 

(66) 
(67) 


dwdfi 


(68) 

(69) 


where  the  Jacobian  is  evaluated  at  the  bifurcation  point.  0  is  a  complex  variable  which 
can  be  written  as  0  =  \aPe  where  (from  equation  62) 

a  =  pPlrCL  +  /?2r«3 

9  =  pPu  +  p2ia2. 

This  leads  to  the  following  algorithm 

1.  calculate  /i,  P,  Q  and  uj 

2.  compute  i?(p,p)  and  i?(p,p) 

3.  solve 


and 


Az0  =  ~^B(  p,p) 

(2 iul  ~  A)z2  =  ^B(p,p) 


for  zq  and  z2. 

4.  compute  B(p,zo),  B(p,z2)  and  C(p,p,p) 

5.  compute  ?'qr(.Bp) 

6.  evaluate  P\  and  P2 

7.  perturb  [i  and  solve 

a  =  pPird  +  /J2ra3 

9  =  PPu  +  P2i.a2. 

until  convergence 

8.  calculate  w,  using  the  expressions  for  wq  and  w2.  setting  e  =  1 


B  Matlab  script  for  calculating  POD  modes 

7.7.  read  m  mode  shapes  row  wise  into  u 
ut=transpose (u) ; 

C=u*ut ; 

C=C/m; 

[a, lambda]  =  eigs(C); 
phi=  ut*a; 

phi=phi/ sqrt (m*lambda) ; 

U  write  out  mode  shapes  from  columns  of  phi 
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C  Mat  lab  script  for  calculating  the  Jacobian  Con¬ 
tributions  from  MUSCL 

for  (1=0; 1<5; 1++)  { 

dl=c2 [1] -cl  [1] ; 
d2=c3 [1] -c2  [1] ; 
d3=c4[l]-c3[l]  ; 

den=0 . 5/ (d2*d2+dl*dl+str2) ; 
lim=(d2*dl+str) * (d2+dl) *den; 

v [1]  =c2 [1] +lim  ; 

vld [0+1*3] =dl*4 . 0*lim*den  -  (d2+dl) *d2*den  -  (d2*dl+str) *den; 
vld [1+1*3] =1.0  -  lim*4 . 0*den* (dl-d2)  +  (d2-dl)*(dl+d2)*den; 
vld [2+1*3]  =-d2*4 . 0*lim*den  +  (d2+dl) *dl*den  +  (d2*dl+str) *den; 

den=0 . 5/ (d2*d2+d3*d3+str2) ; 
lim= (d2*d3+str) * (d2+d3) *den; 

v[l+5]=c3[l]  -  lim; 

vld [0+ ( 1+5) *3] =-d2*4 . 0*lim*den  +  (d2+d3) *d3*den  +  (d2*d3+str) *den; 
vld [1+ (1+5) *3] =1 . 0  +  lim*4 . 0*den* (d2-d3)  -  (d3-d2) * (d3+d2) *den; 

vld [2+(l+5) *3] =d3*4 . 0*lim*den  -  (d2+d3) *d2*den  -  (d2*d3+str) *den; 
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